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Wc present an effective atomic interaction potential for crystalline q;-AI2 03 generated by the program potfit. 
The Wolf direct, pairwisc summation method with spherical truncation is used for electrostatic interactions. 
The polarizability of oxygen atoms is included by use of the Tangney-Scandolo interatomic force field approach. 
The potential is optimized to reproduce the forces, energies and stresses in relaxed and strained configurations 
as well as {0001}, {1010} and {1120} surfaces of AI2O3. Details of the force field generation are given, and 
its validation is demonstrated. We apply the developed potential to investigate crack propagation in a-Al203 
single crystals. 



I. INTRODUCTION 

Aluminum oxide is the most commonly used ceramic 
in technological applications. Due to its insulating prop- 
erties, it is frequently adopted in microelectronic devices, 
i.e., field effect transistors, integrated circuits or super- 
conducting quantum interference devices. Another im- 
portant field of application is the coating of metallic com- 
ponents. Alumina covering aluminum is known to pre- 
vent further oxidation of the metal. Since alumina is 
characterized by a high melting point and also a high de- 
gree of hardness, applications at high temperatures and 
high mechanical demands are possible. Mechanical prop- 
erties such as tensile strengths and fracture processes are 
of high importance regarding these applications. Because 
surfaces and internal interfaces play an important part in 
these technologies, numerical modelling is often focused 
on such systems. 

Due to this technological significance, alumina and 
metal/alumina interfaces were frequently investigated 
experimentally-'^"'^ and by ab-initio methods.''"'' First- 
principles methods are well established for calculating the 
work of separation or the surface energy. However, the 
investigation of dynamic processes such as crack prop- 
agation requires systems with significantly more atoms 
and larger timescales than ab-initio based methods nowa- 
days can deal with. Hence, it is an important challenge 
to develop a suitable interaction force field for classical 
Molecular Dynamics (MD) simulations. 

There exist interaction potential approaches for alu- 
mina, which reach a good accuracy by using angular de- 
pendent terms. Pair potentials are also widely applied 
in simulations of alumina^'^""'^ by reason of usually much 
lower computational demands. Due to the lattice misfit, 
especially the simulation of interface structures has to be 
performed with many atoms. Hence, for the simulation 
of crack propagation, a pair potential is to be favored 
concerning the computational effort. 

Many-body effects can be important for correctly de- 
scribing bond angles and bond-bending vibration fre- 



quencies in oxides. However, for Si02^^"^^ and MgO,^^ it 
was shown that an accurate description can be achieved 
without angular dependent terms, but with polarizable 
oxygen atoms. There, the potential model of Tangney 
and Scandolo^^ (TS) is applied, in which the dipole mo- 
ments of oxygen atoms are determined self-consistently 
from the local electric field of surrounding charges and 
dipoles. 

In Ref. 16, the computational effort in simulations 
scales linearly in the number of particles, which is due 
to the Wolf summation."''' We apply this combination of 
the TS model and the Wolf summation method in the po- 
tential generation as well as in simulations. An effective 
atomic interaction potential is generated with the pro- 
gram potfit. The potential parameters are optimized 
by matching the resulting forces, energies and stresses 
with respective ab-initio values with the force matching 
method. Simulations are performed with the MD code 
IMD.20 

A detailed description of the adopted methods and the 
potential generation is given in Sec. II. MD simulations 
and ab-initio calculations which were performed in or- 
der to validate the obtained force field are presented in 
Sec. III. As a first appplication of the new potential, we 
investigate crack propagation of alumina single crystals. 
To our knowledge, these are the first MD simulations of 
crack propagation in alumina with a potential that takes 
into account the polarizability of oxygen atoms. Results 
of these simulations arc shown in IV. Finally, conclusions 
and an outlook are presented in Sec. V. 



II. FORCE FIELD GENERATION 

A. Tangney-Scandolo force field 

The TS force field is a sum of two contributions: a 
short-range pair potential of Morse-Stretch (MS) form, 
and a long-range part, which describes the electrostatic 
interactions between charges and induced dipoles on the 
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oxygen atoms. The MS interaction between an atom of 
type i and an atom of type j has the form 



MS 



D, 



2exp[f (1-^)] 



(1) 

with r.ij = {fijl, Tij ~ Vj — Vi and the model parameters 
Dij , and rf^ , which have to be determined during the 
force field generation. 

Because the dipole moments depend on the local elec- 
tric field of the surrounding charges and dipoles, a self- 
consistent iterative solution has to be found. In the TS 
approach, a dipole moment pf at position in iteration 
step n consists of an induced part due to an electric field 
E{ri) and a short-range part pf^ due to the short-range 
interactions between charges qi and qj . Following Rowley 
et al.,"^^ this contribution is given by 
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hij and are parameters of the model. Together with 
the induced part, one obtains 



p'l = a,E{r,- {p" 'b=i.jv, {r,b=i.jv) +p: 
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where ai is the polarizability of atom i and E{ri) the 
electric field at position r^, which is determined by the 
dipole moments Pj in the previous iteration step. Tak- 
ing into account the interactions between charges U^"^, 
between dipole moments C/pp and between a charge and 
a dipole J/p'^, the total electrostatic contribution is given 
by 



PP 



and the total interaction is 
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B. Wolf summation 

The electrostatic energies of a condensed system are 
described by functions with r"" dependence, n G 
{1,2,3}. For point charges it is common to ap- 

ply the Ewald method, where the total Coulomb energy 
of a set of N ions, 
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is decomposed into two terms E^"^ and E'^'^ by inserting 
a unity of the form 1 — erfc(Kr) 4- erf (rer) with the error 
function 



erf(Kr) 



dt( 



(8) 



The splitting parameter k controls the distribution of 
energy between the two terms. The short-ranged erfc 
term is summed up directly, while the smooth erf term 
is Fourier transformed and evaluated in reciprocal space. 
This restricts the technique to periodic systems. How- 
ever, the main disadvantage is the scaling of the compu- 
tational effort with the number of particles in the sim- 
ulation box, which increases as 0(Af'^/^),^^ even for the 
optimal choice of k. 

Wolf et al}'' designed a method with linear scaling 
properties 0{N) for Coulomb interactions. By taking 
into account the physical properties of the system, the 
reciprocal-space term E'^ is disregarded. It can be writ- 
ten as 
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where S[k) is the charge structure factor and V the vol- 
ume of the simulation box. The charge structure factor 
is the Fourier transform of the charge-charge autocorre- 
lation function. In condensed matter, the charges screen 
each other; the approach yields proper results in liquids 
and largely also in solids. This means that for small mag- 
nitudes of the wave vectors fc, the charge structure factor 
is also small. In previous studies^^'^'^ it was shown that 
the reciprocal-space term is indeed negligible compared 
to the real-space part, provided that the splitting param- 
eter K is chosen small enough. As however a smaller k 
results in a larger real-space cutoff Tc, the latter has to 
be chosen three or four times the size of a typical short- 
range cutoff radius in metals. In a recent study^^ a cutoff 
radius of only 8 A was sufficient to describe the electro- 
statics of liquid silica and magnesia accurately. 

In addition, a continuous and smooth cutoff of 
the remaining screened Coulomb potential i?'^'^(r,j) = 
qiqj CTk{Krij)r^j^ is adopted at Vc by shifting the poten- 
tial so that it goes to zero smoothly in the first two deriva- 
tives at r = Tc. We use the Wolf method for Coulomb 
and dipolar interactions. For more information about the 
Wolf summation of dipole contributions, the estimation 
of errors attended by this method and a detailed analy- 
sis of the energy conservation in MD simulations see Rcf. 
23. 



C. Generation with potfit 

The programm potfit generates an effective atomic in- 
teraction force field solely from ab-initio reference struc- 
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tures. The potential parameters are optimized by match- 
ing the resulting forces, energies and stresses to according 
first-principles values with the force matching method. 
All reference structures used in this study were prepared 
with the first-principles plane wave code VASP.^^'^^ For 
Nm particles, reference configuration m provides one en- 
ergy ej^, six components of the stress tensor sj^ j {I = 
1, 2, .., 6) and iNm total force cartesian components „ 
(n = 1, 2, iNjn) on A^,„ atoms. The function 

Z = We^e + WsZs + Zf (10) 

is minimized, where 

M 

m— 1 

M 6 

m=l 1=1 

M 3N„, 

m— 1 n— 1 

and Cm, Sm_/ and /m,n are the corresponding values cal- 
culated with the parametrized force field. We and Ws are 
weights to balance the different amount of available data 
for each quantity. They are defined in such a way, that 
each weight can be taken as the pcrcentual amount of 
the concerning force data. In the following we assume M 
reference structures that all consist of the same number 
of particles {N^ = iV), but in principle, potfit can han- 
dle any different number of particles for each reference 
structure. The root mean square (rms) errors, 

AF, = , AFs = and AFf = ./^, 

V iMN V MN ' V MN ' 

(12) 

are first indicators of the quality of the generated force 
field. Their magnitudes are independent of weighting 
factors, number and sizes of reference structures. For 
the minimization of the force field parameters, a com- 
bination of a stochastic simulated annealing algorithm^^ 
and a conjugatc-gradient-likc deterministic algorithm^^ 
is used. For more information about the implementation 
of long-range electrostatic interactions including the TS 
force field with Wolf summation see Ref. 16. 



D. Parametrization 

We first prepared a set of 67 a-alumina crystal struc- 
tures with 360 atoms, in total 24 120 atoms. This refer- 
ence database is composed of three kinds of structures: 
(i) crystals strained up to 20% in [0001], [2110] and [0110] 
directions at zero Kelvin, (ii) structures with free (0001), 
(2110) and (OllO) surfaces at zero Kelvin and (iii) equili- 
brated snapshots taken out of an ab-initio MD trajectory, 
where an ideal crystal is heated from zero Kelvin up to 
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ao 


bAl-O 


CAl-O 


1.122 608 


-0.748 406 


0.026 576 


18.984 286 


-5.571 329 
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Al-Al 




0.000 890 


12.737 442 


5.405 175 


Al-O 




1.000 058 


8.077 778 


1.851 806 


0-0 




0.005 307 


12.081 851 


3.994 815 



TABLE I. Force field parameters, given in IMD units set eV, 
A and amu (hence charges are multiples of the elementary 
charge) . 



2000 Kelvin. It has to be mentioned that no initial ad- 
hoc potential as used in Rcfs. 14 and 16 is required to 
generate a reference database. This is only necessary in 
the case of liquid structures. The present work, how- 
ever, deals with crystalline structures where a database 
can be generated without need of classical MD trajecto- 
ries. The 67 structures are applied as input configura- 
tions for the first-principles plane wave code VASP.^^'^^ 
We used ultrasoft pseudopotcntials^* and the Pcrdew- 
Wang 91^^ generalized gradient approximation (GGA) 
to the exchange-correlation functional. The latter is rec- 
ommended in Ref. 16. A plane wave cutoff of 396 eV was 
used. We adopt a gamma centered k mesh of 2 x 2 x 2 
k points. The weights in potfit were chosen to We — 0.03 

and Ws = 0.28. Setting k = 0.1 A , a cutoff radius of 10 
A was found to be sufficient. The final set of parameters 
is shown in Table I. The rms errors are AF^ = 0.0492 
eV, AFs = 0.0273 eVA-^ and AFf = 0.3507 eVA-^. 



III. VALIDATION 

Firstly, we determine basic properties of crystalline 
alumina such as lattice constants, cohesive energies and 
vibrational properties. Additionally and relevant for frac- 
ture studies, our validation simulations focus on surface 
relaxations, surface energies and stresses of strained con- 
figurations. Finally, we probe the new potential beyond 
its optimization range by investigating two basal twins. 

Lattice constants and cohesive energies are obtained 
by pressure relaxation: Besides energy minimization, the 
pressure tensor of the sample is calculated at each step, 
and the size of the simulation box is changed by a small 
amount in order to lower that pressure. After inserting 
surfaces or interfaces into the relaxed sample, a further 
relaxation is performed which reveals surface or interface 
energies. The ab-initio calculations of these properties 
are performed with VASP, as described in section II D. As 
can be seen from Table II, the lattice constants and the 
cohesive energy of crystalline a-alumina obtained with 
the new potential agree well with our ab-initio calcula- 
tions and previous studies. 

The partial vibrational density of states (VDOS) for 
Gai(E) and Go(E) was obtained by computing the 
Fourier transform of the time-dependent velocity- velocity 
autocorrelation function from a 100 ps MD trajectory 
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E{ii20} 


E{ioio} 


Etwin 


Etwin 




(A) 


(A) 


(cV) 


(J/m^) 


(J/m^) 


(J/m^) 


(J/m2) 


(J/m^) 


New potential 


4.79 


12.97 


31.85 


1.59 


1.65 


1.89 


1.32 


0.57 


Ab-initio 


4.78 


13.05 


32.31 


1.55 


1.83 


1.98 


2.12 


0.76 


Lit. (experiment) 


4.76-'" 


is.oo-"" 


31.8^' 












Lit. (ab-initio) 






33.0^^ 


1.85''^ 


2.44^^ 


2.39^^ 


1.99" 


0.73* 










2.03^* 


2.23^* 


2.50^^ 







1.76' 



TABLE II. Lattice constants, cohesive energies and interface energies of basal twins obtained with the new potential compared 
to ab-initio results and literature data. 



with the software paekage nMoldyn.^ The generalized 
VDOS G(E) is then calculated by 

G{E)^ ^ ^G,{E) (13) 

with the scattering cross section cr^ and atomic mass 
of atom fi. Fig. 1 shows the partial and generalized 
VDOS. The key features of the curves obtained with the 
new potential and an ab-initio study from Ref. 38 coin- 
cide. For the partial VDOS of aluminum, both studies 
show a broad band between 41 and 83 meV. The ab- 
initio results show less states in the low-energy band be- 
tween 15 and 30 meV. There are two sharp peaks between 
87 and 100 mcV in the ab-inito curve, whereas the new 
potential shows one broader peak with a shoulder, that 
indicated the second peak. The curves for the partial 
VDOS of oxygen are in good agreement. Both simulation 
and ab-initio calculation show the main band of states 
between 12 and 85 mcV with similar curve progression 
and three local maxima at around 35, 47 and 62 meV. 
These characteristics are also reflected in the generalized 
VDOS. In addition, the generalized VDOS obtained with 
neutron scattering^ is depicted, which shows the same 
main characteristics of the curve on a qualitative level. 

The literature surface cnergics'^^"'^^ given in Table II 
differ among each other which origins from the different 
methods used in the ab-initio approaches. Our calcula- 
tion of the (0001) surface energy agrees with the value 
obtained in Ref. 36, where also the VASP code with the 
same pseudopotential and exchange-correlation approxi- 
mation was used. For all investigated surfaces, the ener- 
gies obtained with the new potential and with ab-initio 
methods agree. Both results reveal that the (0001) sur- 
face has the lowest surface energy and the (0110) surface 
the highest one. The studies of Ref. 33-35 yield slightly 
higher energies for all investigated surfaces. 

The surface structure after relaxation is shown in 
Fig. 2. It reveals that the Al atoms of the outermost 
Al-layer are moved closer to the outermost 0-layer at 
the (0001) Al-terminated surface, which is known to be 
the most stable (0001) surface termination. The atomic 
adjustment perpendicular to the surface obtained with 
the new potential agrees very well with our ab-initio re- 
sults. The distance of an Al- and an 0-laycr is 0.83 A, 
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FIG. 1. VDOS calculated with the new potential compared 
with an ab-initio study^* and experiment.^ a) Partial VDOS 
for aluminum atoms, b) partial VDOS for oxygen atoms, c) 
generalized VDOS. 



whereas this value decreases to 0.15 A (MD) respectively 
0.14 A (ab-initio) at the surface. A relaxation of the oxy- 
gen atoms can be seen at the (0110) surface. Both MD 
and ab-initio study show that the three oxygen atoms 
per unit cell - which are initially in a row along the di- 
rection orthogonal to the plane of Fig. 2 - relax to dif- 
ferent distances from the initial surface. Furthermore, a 
relaxation towards the neighboring Al-atoms in the first 
layer occurs. The results obtained with the new poten- 
tial again coincide with the ab-initio calculation. One 
difference, however, can be seen at the second layer: Ev- 
ery second Al atom in each row orthogonal to the figure 
plane is slightly displaced towards the surface in the ab- 
initio relaxation. This effect is not observed in the MD 
relaxation simulation. Neither the MD nor the ab-inito 
relaxation study of the (2110) surface yield significant 
atomic movements. 

As described in section II D. strained configurations 
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FIG. 2. Surface relaxations for the new potential compared 
with ab-imtio calculations. Red: oxygen. Blue: aluminum, 
a) (OOOl)-surface, ab-imto, b) (OOOl)-surface, MD, c) (0110)- 
surface, ab-inito, d) (OllQ)-surface, MD. 

with strain directions [0001], [2110] and [OllO] were 
added in the reference database for the potential opti- 
mization approach. To clarify, whether the new poten- 
tial can reproduce the stresses of strained configurations, 
these stresses are calculated using MD with the new po- 
tential. 

Fig. 3 shows a comparison of stresses obtained in sim- 
ulations to the stresses of the underlying ab-initio ref- 
erence configurations, each with 360 atoms. They are 
strained up to 15% in [0001], [2110] and [0110] direction 
respectively. The difference between MD and ab-initio 
results is small at lower strains. With increasing strain, 
the difference increases up to about 10 GPa at 15% strain. 
The new potential underestimates the stress in all cases. 
However, the directions, in which the stress increases, 
can be reproduced correctly: The highest stresses are ob- 
served for strains in [0001] direction, the lowest stresses 
are found for configurations strained in [0110] direction. 

Finally, we probe the potential by simulating basal 
twins. These are systems which were not included in 
the optimization process. This demonstrates the appli- 
cability of the new force field beyond the range for which 
it was optimized. It is well known that the basal twin 
is a frequently observed AI2O3 interface. We consider 
two types of basal twins which differ by their interface 
structure: The mirror twin consists of two single crys- 
tal elementary cells with mirror symmetry. The mirror 
plane is a (0001) oriented oxygen layer. Based on the 
mirror twin, one can obtain the rotation twin by shifting 
one of the single crystal elementary cells in the [0110] 
direction. For details of the interface structures see Ref. 
4. As can be seen from Table II, the ab-inito calculated 
interface energies of the basal twins agree well with the 
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FIG. 3. Stresses of strained configurations with [0001] (green), 
[2110] (red) and [0110] (blue) strain directions against the 
stress tensor component in direction of strain. The curves 
obtained with the new potential are marked with squares, 
those from ab-initio calculations with circles. 



results from Ref. 4. The MD simulations with the new 
potential underestimate both energy values. But the po- 
tential is able to qualitatively reproduce the relation of 
the two interface energies: The interface energy of the 
mirror twin (E"4in) about two by three times the one 
of the rotation twin (E^^jj^); that is why, the latter one 
is energetically favorable. 



IV. CRACK PROPAGATION 

A. Simulation conditions 

To investigate crack propagation in mode I at constant 
energy release rate, we prepared configurations with di- 
mensions {bx,by,bz) of 21 X 3 X 13 nm'^ which contain 
about 80 000 atoms. Lateral planes of these cuboids are 
the (0001), (0110) and (2110) crystallographic planes. 
An elliptical initial crack of 5 nm length in x-direction 
is inserted on one side of the samples by moving atoms 
in z-direction. The opening of the crack is calculated 
with the Griffith criterion: In front of the crack tip the 
sample is linearly strained until the elastic energy due to 
strain is equal to the Griffith energy Go which is twice 
the surface energy of the crack plane. This criterion was 
fuUfiUed at about 4% strain. These structures are then 
relaxed to obtain the displacement field of a stable crack. 
Periodic boundary conditions are applied in the direc- 
tion along the crack front, whereas fixed displacement 
boundary conditions arc applied in the other directions. 
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FIG. 4. Crack propagation in Q-AI2O3 (Al-atoms blue, O- 
atoms red), a) Initial crack in (2110) plane, b) initial crack 
in (0001) plane, c) initial crack in (0110) plane. 



Initial configurations for crack propagation simulations 
with different energy release rates are obtained by a lin- 
ear scaling of this displacement field. During the crack 
propagation simulations, the NVE (constant number of 
particles, volume and energy) ensemble and a starting 
temperature of K are applied. A timcstcp of 1 fs is 
used for all simulations. 



not the preferred fracture planes. 



C. Cracks in a {0001} plane 

The {0001} planes are the closest packed planes of a- 
AI2O3; therefore cleavage of these planes is improbable. 
Our simulations confirm that cracks do not propagate in 
the (0001) plane in which the initial crack was inserted 
with crack propagation directions [2110] or [0110]. At en- 
ergy release rates below 1.7 Gq, a damaged region in front 
of the crack tip is generated with atomic disorder simi- 
lar to the case of the (2110) plane. However, the crack 
tip does not move. More interesting is the (0001) [0110] 
oriented initial crack at higher energy release rates: As 
can be seen from Fig. 4 b, the crack propagates in a 
{IOI2} cleavage plane. Due to the boundary conditions, 
the crack surfaces cannot separate completely, but a row 
of oxygen atoms stays in the crack path. In contrast 
to the cracks in a {1120} plane, there is no disorder 
along the crack surfaces. The observed crack propaga- 
tion in a {1012} cleavage plane agrees well with electron 
microscopy investigations^ which revealed that fracture 
surfaces of the {1012} cleavage planes are frequently ob- 
served in alumina. 



D. Cracks in a {1010} plane 



B. Cracks in a {1120} plane 

Initial cracks were inserted in a (2110) plane with crack 
propagation directions [0001] or [0110]. In both cases, 
the initial crack is stationary at energy release rates up 
to 1.5 Go . Higher energy release rates lead to crack prop- 
agation in the initial (2110) plane. It is known that the 
discrete character of interatomic bonds which reveals it- 
self in the so called lattice trapping effect^'' can retard 
crack propagation to energy release rates beyond the crit- 
ical energy release rate determined by the Griffith crite- 
rion. As can be seen from Fig. 4 a, yet at high energy 
release rates bond breaking between the two crack lips is 
not continuous. Chains of atoms bridging the crack lips 
stay intact during the simulation. Furthermore, consid- 
erable disorder is observed at the resulting crack surfaces. 
In [0110] propagation direction, these effects are slightly 
more pronounced than in [0001] propagation direction. 
It cannot be ruled out that these very strong bonds arc 
caused by inaccuracies of the force field. However, the 
observation of intact bonds bridging the crack lips agree 
with previous simulations of alumina, where it was stated 
that these atomic chains were observed in experiments as 
well.^'^ Our results show that cracks are able to propa- 
gate in (2110) planes in both [0001] and [0110] directions. 
This is in accordance with experiments,^ where it turned 
out that cracks occur in {1120} planes although these are 



Cracks initially inserted in a (0110) plane in [2110] di- 
rection propagate in this orientation starting at energy 
release rates of 2.2 Gq. In the case of a [0001] crack propa- 
gation direction (Fig. 4 c), the movement of the crack tip 
was observed at energy release rates above 1.9 Gq. Simi- 
lar to the case of the (0001) [0110] orientation, the crack 
changes its propagation plane. As can be seen in Fig. 4 
c, the crack moves partially in the initial (0110) plane, 
but also partially in a {1012} cleavage plane. As in the 
case of cracks in {1120} planes, the initially (0il0)[2110] 
or (0110) [0001] oriented cracks generate disorder at the 
crack surfaces during crack propagation; some atomic 
bonds stay intact across the crack opening. 



E. Orientation of electric dipole moments 

The orientation of electric dipole moments in crack 
propagation simulations can differ significantly from the 
orientations in other configurations of the same material. 
There are two main effects which influence the align- 
ment of dipoles: Charged surfaces and piezoelectricity. 
In order to investigate these effects in a-Al2 03, dipoles 
of the crack propagation simulations are visualized with 
the program MegaMol.^"'^ Since effects from the boundary 
conditions (which involve unnatural surfaces in some di- 
rections) should be excluded, we performed a crack prop- 
agation simulation with periodic boundary conditions in 



FIG. 5. Electric dipole moments of oxygen atoms in crack propagation simulation. Each oxygen atom is visualized by an arrow 
with the direction of the dipole moment; the lengths are normalized. The colour corresponds to the orientation: blue - down, 
yellow - up, green - left, red - right. 



all directions. In this case, a symmetrical crack was in- 
serted in the middle of the sample. As long as the two 
crack tips are far away from each other (which is true up 
to about 5 nm distance), we observe the same result as 
with the boundary conditions described above. Hence, 
possible boundary condition effects can be excluded. 

Charged crack surfaces occur as a function of the crack 
propagation planes. As shown in Fig. 4, the initial sur- 
faces of cracks in a {1120} plane are oxygen terminated. 
Therefore, it can be expected that the dipoles of the oxy- 
gen atoms at each crack surface are aligned in the same 
direction and the alignment on the other crack surface is 
the other way round. This adjustment is reproduced in 
our simulation, as can be seen in Fig. 5. It is significantly 
more pronounced at the surfaces of the initial crack (left 
part of the crack in Fig. 5), but it can be seen also along 
the crack surfaces created by the propagation (right part 
of the crack). In the latter case, the disorder along the 
crack surfaces and the intact atomic chains across the 
crack opening reduce the alignment. 

Collective orientation mechanisms of electric dipole 
moments due to strain can only be observed in crystals, 
where the inversion symmetry is broken. Therefore, in 
bulk a-Al203 usually no dipole alignment is observed. 
With the new force field, we simulated a-Al203 strained 
in [2110] direction up to 15%. No collective alignment of 
dipole moments occurcd in this simulation. In the crack 



propagation simulation, however, the inversion symmetry 
is broken by the crack. As can bee seen from Fig. 5, an 
orientation of the electric dipole moments due to strain 
is observed. Dipoles are aligned in the region in front 
of the crack tip which corresponds to the strained part 
of the configuration. On the contrary, in the unstrained 
regions below and above the crack surfaces the dipoles 
show no preferred orientation. 



V. CONCLUSIONS 

We presented the generation and validation of an ef- 
fective ab-initio based polarizable force field for MD sim- 
ulations of a-Al2 03. It was shown that it can be used to 
simulate mechanical and vibrational properties of crys- 
talline alumina with high accuracy. Furthermore, the 
potential is designed for the simulation of surfaces, no- 
tably for crack propagation. Due to the Wolf summa- 
tion, we achieve linear scaling properties in the number 
of particles, which makes it possible to investigate typi- 
cal required system sizes and time scales for crack simu- 
lations in resonable real time. The polarizability of oxy- 
gen atoms is taken into account by use of the Tangney- 
Scandolo interatomic force field approach. To our 
knowledge, the present work is the first study of the be- 
havior of electric dipole moments in MD simulations on 
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crack propagation. 

As a first application beyond the study of basic crys- 
talline properties, we performed crack propagation sim- 
ulations with initial cracks in different crystallographic 
planes of q:-A12 03. Wc have shown that cracks usually 
propagate in the initial plane, but no crack propagation 
occurs in the close packed {0001} planes. It was also 
observed that cracks tend to deflect to a {10l2} cleav- 
age plane in the case of an initial crack plane which is 
unfavourable regarding crack propagation. The simula- 
tion result, that {1012} cleavage planes are favourable 
regarding crack propagation, agrees well with electron 
microscopy investigations^ of cracks in a-Al2 03. Addi- 
tionally, it was shown that the new force field allows to 
investigate electric dipole orientations in various strained 
structures. Dipole alignment due to strain in a-Al2 03 
was found in the performed crack propagation simula- 
tions. 

In the future, the new force field can be applied for 
large scale simulations of alumina systems, particularly 
with regard to surfaces and interfaces, e.g. polycrystals, 
or multi-layer composites. The new force field allows 
further increasing of system sizes to several millions of 
atoms at time scales up to 1 ns. 
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